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An algorithm for efficiently incorporating the effects of mesh sensitivities in a computa- 
tional design framework is introduced. The method is based on an adjoint approach and 
eliminates the need for explicit linearizations of the mesh movement scheme with respect to 
the geometric parameterization variables, an expense that has hindered practical large-scale 
design optimization using discrete adjoint methods. The effects of the mesh sensitivities can 
be accounted for through the solution of an adjoint problem equivalent in cost to a single 
mesh movement computation, followed by an explicit matrix-vector product scaling with the 
number of design variables and the resolution of the parameterized surface grid. The accu- 
racy of the implementation is established and dramatic computational savings obtained using 
the new approach are demonstrated using several test cases. Sample design optimizations are 
also shown. 


Nomenclature 

D = vector of design variables 

/ = cost function 

K = linear elasticity coefficient matrix 

L = Lagrangian function 

L/D = Lift-to-drag ratio 

Q = vector of dependent variables 

R = discretized residual vector 

X = computational mesh 

A = vector of adjoint variables 

I. Introduction 

In recent years a concerted effort has been made to bring higher fidelity, physics-based computational fluid dy- 
namics (CFD) simulations into the aircraft design process. Such tools have typically been used to perform validations 
of designs derived through the use of lower fidelity methods, or in some cases, used in heuristic design methods that 
rely heavily on experience. However, these advanced CFD tools are now routinely targeted as primary components of 
automated optimization frameworks. 

In the field of gradient-based design, one challenge in utilizing solvers based on the Euler or Reynolds-averaged 
Navier-Stokes equations has been an efficient and accurate method for computing the sensitivity derivatives required 
by many optimization algorithms. Approaches such as finite differencing, 1 * " 3 direct differentiation, 4 ' 10 and the com- 
plex variable method 1116 can be used for calculating these derivatives; however, their cost scales directly with the 
number of design variables. For typical aerodynamic design problems where this value may be on the order of tens to 
hundreds, this limitation precludes the use of such methods. 

To alleviate the computational burden associated with problems containing many design variables, recent work has 
focused on the use of adjoint methods. 17 " 35 Adjoint methods may be implemented in either a continuous or discrete 
context, depending on the order in which the differentiation and discretization operations are performed. The relative 
merits of the two approaches are the subject of much debate in the literature; one such difference is the focus of the 
current work. 

Both the continuous and discrete adjoint approaches introduce an auxiliary, or adjoint, variable that is determined 
through the solution of an additional linear system of equations. Using the same solution procedure as is used for the 
governing equations, asymptotic convergence rates of both systems are similar, as the eigenvalues of the two systems 
are the same. Moreover, for the discrete adjoint variant, if the solution algorithm itself is constructed in a manner 
which is discretely adjoint to the baseline scheme, the asymptotic convergence rates are guaranteed to be identi- 
cal 19,20,31 
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The primary advantage of the adjoint approach is that the expense does not scale with the number of design vari- 
ables. Indeed, the solution of the adjoint system itself is independent of the number of design variables. In practice, 
however, the discrete approach dictates that the effects of the mesh be accounted for in the sensitivity derivatives, 
whereas these terms do not appear in the derivation of the continuous operators. While perhaps counterintuitive, this 
difference is actually viewed as one of the primary advantages of the discrete approach where, by definition, the sen- 
sitivity analysis is discretely consistent with the flow simulation. If these grid-related contributions are not included, 
it has been shown that the resulting sensitivities can be extremely inaccurate and even of the incorrect sign. 29,30 The 
continuous approach should converge to the discrete result in the case of a sufficiently refined grid, but this condition 
is seldom, if ever, met in current practice for realistic three-dimensional complex geometry computations. 

Unfortunately, the expense of computing these grid-related terms has hindered large-scale application of the dis- 
crete approach. For a single design variable, the use of a direct mode of differentiation to obtain the derivative of a 
mesh movement scheme based on linear elasticity was shown in Ref. 30 to cost as much as 30% of a flow simulation. 
Clearly, the expense associated with performing this operation for several dozen design variables will dominate the 
overall computation. 

In the current work, an adjoint formulation is introduced that obviates the need for an explicit computation of the 
grid sensitivities. Such an approach has previously been used in the automatic differentiation (AD) community; 36 
however, no explicit formulations for handcoding such a technique have been found in the literature. The scheme pre- 
sented in the current work significantly reduces the expense associated with the mesh linearizations. Demonstrations 
of the method show that a rigorous discrete sensitivity analysis for problems based on the Navier-Stokes equations 
may ultimately be performed at a cost comparable to that of the analysis problem. 

II. Mesh Sensitivities via Forward Mode Differentiation 

The discrete adjoint technique for sensitivity analysis can be derived in several ways. Here, the approach taken in 
Ref. 28 is used. Consider the vector of discretized residual equations R for the Euler or Navier-Stokes equations as a 
function of the design variables D , computational mesh X , and flowfield variables Q . Given a steady-state solution 
of the form R(D, Q, X) = 0 , a Lagrangian function L can be defined as 

L(D, Q , X, A) = f(D, 0, X) + A T R(D, Q , X) (1) 

where f(D , Q, X) represents a cost function to be minimized and A represents a vector of Lagrange multipliers, or 
adjoint variables. Differentiating this expression with respect to D yields the following: 
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Since the vector of adjoint variables is essentially arbitrary, the coefficient multiplying [dQ/dD] can be eliminated 
using the following equation: 


= Jf 
La0j dg 


(3) 


Eq. 3 represents the discrete adjoint equation for the flow simulation. The solution of this linear system of equations 
for three-dimensional turbulent flows on unstructured grids has been demonstrated previously in Refs. 28-31. Once 
the solution for A has been formed, the remaining terms in Eq. 2 can be evaluated to give the desired sensitivity vec- 
tor: 
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The [3X/3DJ terms in Eq. 4 represent the mesh sensitivities. In Ref. 30, a mesh movement strategy based on the 
equations of linear elasticity is described. Although generally not as expensive as a flowfield computation, the cost 
associated with solving these equations is not trivial. If the system is posed as 

KX = X surface , (5) 

then the mesh sensitivities may be computed from the following: 
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Note that the solution of this linear system is equivalent in cost to that of the mesh movement (Eq. 5), and must be ob- 
tained once for each design variable in D . 

III. Adjoint Approach for Eliminating Mesh Sensitivities 

Now reconsider Eq. 1, where a subscript/has been appended to A to indicate the adjoint variable associated with 
the flow equations, and an additional adjoint variable A multiplying the residual of the grid movement problem has 
been introduced: 


L(D, Q , X, A f , A,,) = f(D, Q , X) + aJ R(D, Q , X) + A g (KX - X surface ) (7) 

Linearizing with respect to D as before yields 


dL 

dD 


5/ 

dD 


dR 

dD 


-J 


A/ + 


dQ 

dD 


-J 


dR 


3 / 

dQ ' IdQj 


-J 


dX 


-J 


dD 


3/ 

dX 


dR 

dX 


-J 


A/+a; 


K \- A 


raz- 

dD 


Asurface 


(8) 


j 

As before, the coefficient multiplying [dQ/dD\ can be eliminated by satisfying the adjoint equation for the flow 
simulation, Eq. 3. In a similar fashion, the term multiplying [ dX/dD ] can also be eliminated by satisfying a second 
adjoint problem: 


K T A„ 




( 9 ) 


With the solution of Eqs. 3 and 9, the final form of the sensitivity vector becomes 
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With the formulation outlined above, a single solution of the linear system given by Eq. 9 is required for each func- 
tion /. If / is some quantity composed of aerodynamic coefficients such as lift, drag, or moments, several observa- 
tions can also be made about Eq. 10. For design problems in which the shape is held constant and only global param- 
eters such as the angle of attack or Mach number are allowed to vary, Eq. 9 need not be evaluated, and the third term 
in Eq. 10 is identically zero. Conversely, for problems involving solely geometric parameterization variables, the first 
and second terms in Eq. 10 are identically zero, as there is no explicit dependence of / or R on D . In addition, the 
term A g [dX/dD] surface ' s ver Y cheap to compute and only requires an explicit inner product dimensioned by the 
size of the surface mesh for each design variable. 

Implementation 

The software developed in Refs. 28-30 has been extended to include the formulation outlined above. In the previ- 
ous implementation, the matrix-vector product [dR /dX\ A ( has been stored to enable rapid computation of each 
inner product with dX/dD , as these residual linearizations are constant for all shape parameters in I) . These me- 
chanics are now used to construct the right-hand side of Eq. 9 in the current work. 

The parameterization schemes employed here are described in Refs. 37 and 38, and rely on a free-form deforma- 
tion technique to provide a compact set of design variables for a wide range of configurations. Given the current vec- 
tor of design variables D , the methods are used to determine the current location of the surface grid points as well as 
their analytic derivatives with respect to D , [dX/dD] sur ^ ace . The schemes are very inexpensive and can be used to 
consistently parameterize families of computational meshes suitable for multiple disciplines. 

The matrix K is formed using the method of Ref. 30, and is merely transposed once all of the contributions have 
been included. The system given by Eq. 9 is then solved using the Generalized Minimal Residual (GMRES) 
algorithm 39 implemented in Ref. 30. Since the eigenvalues of K remain unchanged by the transpose operation, the 
solution for A g converges similarly to that of the mesh movement scheme. In the current implementation, adjoint so- 
lutions for multiple functions / may be computed simultaneously as outlined in Ref. 3 1 by storing multiple right- 
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Figure 1 ONERA M6 grid used for evaluating 
linearizations. 



Figure 2 Convergence rates for direct and 
adjoint modes. 


hand sides for Eq. 9. Once the solution for A„ has been computed, an explicit inner product with the surface mesh 
sensitivities for each design variable yields the final sensitivity vector. 

IV. Consistency of Linearization 

To verify the accuracy of the implementation, derivatives of the lift and drag coefficient are computed for fully tur- 
bulent flow over the ONERA M6 wing geometry 40 shown in Fig. 1. The grid for this case contains 16,391 nodes and 
90,892 tetrahedra. The freestream Mach number is 0.84, the angle of attack is 3.06 degrees, and the Reynolds number 
is 1 million based on the mean aerodynamic chord. The computations have been performed in a domain-decomposed 
environment on 12 Pentium IV processors using three different approaches as outlined in Table 1. All equation sets 
have been converged to machine accuracy. 

The first linearization method is a direct mode of differentiation using complex variables. This approach was origi- 
nally suggested in Refs. 12 and 13, and was first applied to a Navier-Stokes solver in Ref. 11. The primary advantage 
of this approach is that true second-order accuracy may be obtained by selecting step sizes without concern for sub- 
tractive cancellation error typically present in real-valued divided differences. Through the use of an automated 
scripting procedure as outlined in Ref. 41, this capability can be immediately recovered at any time for the baseline 
code. For computations using this method, the complex step size has been chosen to be 1x10 . The second tech- 

nique used to verify the linearizations relies on the handcoded implementation 29,31 of the discrete adjoint system 
given by Eq. 3 for the flow equations and the handcoded direct differentiation of the mesh terms 30 as formulated in 


Table 1 Schemes used to obtain sensitivities. 


Method 

Flowfield Linearization 

Mesh Linearization 

1 

Direct differentiation with 
automated complex variables 

Direct differentiation with 
automated complex variables 

2 

Handcoded discrete adjoint 
(Eq. 3) 

Handcoded direct differentiation 
(Eqs. 4 and 6) 

3 

Handcoded discrete adjoint 
(Eq. 3) 

Handcoded discrete adjoint 
(Eqs. 9 and 10) 
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Table 2 Comparison of sensitivity derivatives for lift and drag coefficients using various approaches. 


Objective 

Method 


Design Variable 


Function 

Thickness 

Shear 

Camber 

Twist 



1 

-0.584383430968430 

-0.073891855284066 

1.843734584180741 

-0.022010251214990 

c L 

2 

-0.584383430967291 

-0.073891855283698 

1.843734584180810 

-0.022010251215005 


3 

-0.584383430968115 

-0.073891855283921 

1.843734584180955 

-0.022010251214989 


1 

0.058894900355748 

-0.006835640271421 

0.064393773359690 

-0.001817294278046 

c D 

2 

0.058894900355806 

-0.006835640271405 

0.064393773359692 

-0.001817294278046 

3 

0.058894900355780 

-0.006835640271392 

0.064393773359720 

-0.001817294278046 


Eqs. 4 and 6. The third approach is similar to that of Method 2; however, the present adjoint formulation given by 
Eqs. 9 and 10 is used to evaluate the mesh terms. 

Sensitivity derivatives for the lift and drag coefficients using several shape parameters located at the midspan of the 
wing are shown in Table 2. All three methods are in excellent agreement, with the only discrepancies apparent in the 
twelfth decimal place or better; note that even at machine precision, the last several digits are often still fluctuating. A 
comparison of the convergence rates for the systems given by Eqs. 6 and 9 is shown in Fig. 2. The systems converge 
in a similar fashion. 


V. Large Scale Test Cases 

Several large scale test cases are used to evaluate the impact of the new scheme. For the cases shown, relative mea- 
sures in performance are used to establish benefits of the current approach. The meshes used here are not necessarily 
optimal for resolving the flowfields; the intent is solely to demonstrate algorithmic performance on a given discrete 
model. For the same reason, vehicle performance benefits obtained here through optimization are not claimed to be 
representative of what may be achieved for grid-resolved computations. Each of the grids shown has been obtained 
using the methods of Refs. 42 and 43. Note that although the cases described here involve compressible flows, an in- 
compressible implementation 29,44 based on the method of artificial compressibility 45 is also maintained. 

To show the impact of the current formulation. Table 3 lists timings for the major software components necessary 
for aerodynamic shape optimization. Operations associated with the surface grid evaluation and sensitivities are omit- 
ted; these costs are negligible in comparison with the other components. The computations for each configuration 
have been performed on different hardware with a range of compilers; timings for a given test case should only be 
compared relative to each other and not to other test cases. For each of the computations, the residuals for Eqs. 5 and 


Table 3 Timings for analysis and sensitivity analysis components. 




Wallclock Time, mins 


Test Case 

Analysis 


Sensitivity Analysis 

(Number of 
Design Variables) 

Mesh 

Movement 

Flowfield 

Solution 

Flowfield 
Adjoint Solution 

Mesh Adjoint 
Solution and 
Sensitivities 

Inviscid Stowed 
Morphing Vehicle 
(175) 

4 

18 

42 

4 

Turbulent Cruise 
Morphing Vehicle 
(150) 

3 

45 

70 

4 

Turbulent 

Wing-Body 

(172) 

8 

228 

184 

10 


5 

American Institute of Aeronautics and Astronautics 






9 have been reduced by ten orders of magnitude. All of the values shown in the table use an objective function based 
on the drag coefficient. The data shown in the table for each test case is discussed in more detail below. 

The time required for a flowfield adjoint solution is discussed at length in Ref. 31. The relative cost of these com- 
putations can vary widely depending on the coupling used for the turbulence closure equation. The storage strategy 
used for the exact linearization of the residual required for the adjoint solution can also have a large impact on the ef- 
ficiency of the adjoint solver. Here, the nearest neighbor terms are stored, while the larger stencil contributions are re- 
computed at each time step in order to save memory. If sufficient memory is available, the entire linearization may be 
stored and the relative cost of the adjoint solution for the flowfield can be drastically reduced by computing an ex- 
plicit matrix-vector product for the residual of Eq. 3. 

In order to save computational time during optimizations, the flow solver and its associated adjoint solver are typi- 
cally restarted from previous solutions, so that, in general, these solutions require less time as the optimization con- 
verges. Therefore the data indicated in Table 3 cannot necessarily be used to directly infer the cost of an entire opti- 
mization. However, the mesh linearizations are recomputed for every sensitivity analysis so the savings of the current 
adjoint formulation over the direct approach scale directly with the number of sensitivity analyses required over the 
course of an optimization. This benefit will be demonstrated in the following sections. 

Morphing Aircraft Configuration 

The first two test cases are based on an unmanned combat air vehicle being pursued by Lockheed-Martin and 
DARPA and is described in Ref. 46. The aircraft is designed to cruise in a conventional configuration and to fold its 
wings into a stowed position for the dash/attack portion of its mission. 

The first test case is based on the stowed configuration and is computed using the Euler equations. The grid used 
for this case is shown in Fig. 3, and contains 430,732 nodes and 2,543,772 tetrahedra. The freestream Mach number 
is 0.80 and the angle of attack is 2 degrees. For this test, the geometry has been parameterized using the package out- 
lined in Ref. 38 to create a set of 175 design variables describing the wing and fuselage. 

Timing results obtained on 16 Pentium IV processors for the various components needed for optimization are 
shown in Table 3. Here, the flow solver and its adjoint counterpart have been run 300 time steps. The analysis por- 
tion of the computation, consisting of a mesh movement and a flowfield solution, requires approximately 22 minutes 
of wallclock time. With a single solution of Eq. 9, the sensitivity analysis for all 175 design variables takes 46 min- 
utes. Note that the previous method based on direct differentiation would require 175 solutions to Eq. 6, or nearly 12 
hours more wallclock time. This savings is the equivalent of 39 additional flowfield solutions. 

To demonstrate an optimization for the current test case, an objective function based on drag has been used in con- 
junction with an explicit lift constraint. In this manner, separate adjoint solutions are required to evaluate the sensitiv- 
ities of both the objective and the constraint. A subset of 63 of the 175 total shape variables are free to change, subject 
to side constraints that prohibit thinning of the geometry. The optimization package described in Ref. 47 is used to 
drive the design problem. The results of the optimization are shown in Fig. 4, where the drag coefficient has been re- 
duced by 52% and the lift coefficient satisfies the requested constraint value. For this case, the optimizer requested 26 



Figure 3 Surface grid for morphing aircraft in dash 
configuration. 



O 

fe 

CD 

O 

o 

o> 


Figure 4 Lift and drag coefficients during 
constrained optimization of dash configuration. 
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Figure 5 Surface grid for morphing aircraft in cruise 
configuration. 



Figure 6 Lift-to-drag ratio during unconstrained 
optimization of cruise configuration. 


sensitivity analyses during the course of the optimization. Based on the timings above, the current adjoint formulation 
results in a savings of roughly 109 hours of wallclock time — or the equivalent of 364 flowfield solutions — in com- 
puting the mesh linearizations for the 63 variables. 

The next case is a turbulent flow computation for the baseline cruise configuration as shown in Fig. 5. The grid 
contains 534,525 nodes and 3,159,677 tetrahedra. Here, the freestream Mach number is 0.80, the angle of attack is 2 
degrees, and the Reynolds number is 5 million based on the mean aerodynamic chord. This case has been run on 64 
Opteron processors, and results for the component timings are again shown in Table 3. For these computations, the 
flow solver and associated adjoint solver have been run for 800 time steps using a loosely coupled formulation for the 
turbulence equation. The mesh linearizations for all 150 design variables have been computed in approximately four 
minutes. With the direct approach, equivalent results would require 150 solutions to Eq. 6 at an expense of roughly 
450 minutes or the equivalent of 10 flowfield solutions. 

For this test case, the optimization goal is to perform an unconstrained maximization of the lift-to-drag ratio (L/D). 
The optimization package used for this computation is described in Ref. 48. Results for the L/D maximization are 
shown in Fig. 6, where a subset of 54 of the 150 shape design variables are allowed to change, subject to side con- 
straints that the geometry does not decrease in thickness. The final value of L/D is approximately 24% higher than its 
initial value. Here, the optimizer required seven sensitivity analyses. The present adjoint implementation for the 
mesh-related terms has yielded an estimated savings of approximately 19 hours of wallclock time, or the cost of 25 
flowfield solutions, over the direct approach. 

Slotted Wing-Body Configuration 

The final example is fully turbulent flow over a transport wing-body configuration with a wing incorporating an 
advanced slotted airfoil concept. The mesh for this case is shown in Fig. 7 and contains 1,326,150 nodes and 
7,744,304 tetrahedra. The freestream Mach number is 0.87, the angle of attack is 2.7 degrees, and the Reynolds num- 
ber is 3 million based on the mean aerodynamic chord. The computations shown in Table 3 have been performed 
using 64 Opteron processors. Here, the flow solver and its adjoint counterpart have each been run 1000 time steps, 
using the tightly coupled implementation of Ref. 31 for the turbulence equation. For this reason, the relative expense 
of the flowfield adjoint solution is less than that of the previous examples, since the complete Jacobian of the residual 
operator is required at every time step during the flow solution. The surface grids for this computation have been pa- 
rameterized using the method of Ref. 37 to create a set of 172 design variables describing the main wing and flap ge- 
ometries. For this case, the analysis problem takes 236 minutes of wallclock time versus 194 minutes for the sensitiv- 
ity analysis. The present adjoint formulation for the mesh terms yields a savings of approximately 23 hours in 
computing the sensitivities for all 172 variables. This is the equivalent of six additional flowfield solutions. 

An unconstrained maximization of L/D is performed for the current configuration, using a subset of 41 of the 
shape parameters in addition to the angle of attack. As in the previous examples, the geometry is constrained to main- 
tain its original thickness. The net result of the optimization procedure is a 32% increase in L/D as shown in Fig. 8. 
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Figure 7 Surface grid for slotted wing-body 
configuration. 



Figure 8 Lift-to-drag ratio during unconstrained 
optimization of wing-body configuration. 


For this case, the optimizer required 26 sensitivity analyses during the course of the computation. The previous for- 
mulation for the grid sensitivities would have required an estimated additional 142 hours of wallclock time to per- 
form the same optimization. This savings is the equivalent of 37 flowfield solutions. 

VI. Summary and Conclusions 

An adjoint formulation to efficiently account for mesh sensitivities has been developed and implemented. In com- 
bination with an adjoint procedure for the aerodynamic flowfield, the new approach effectively removes the depen- 
dence on the size of the design space for a rigorous discrete sensitivity analysis. Unlike previous methods that re- 
quired explicit computation of grid sensitivity terms for each design variable, the current approach accounts for these 
terms through the solution of a single adjoint problem, equivalent in cost to the mesh movement scheme being used. 
The method has been implemented within a three-dimensional unstructured grid framework, and the resulting sensi- 
tivity derivatives have been shown to be in excellent agreement with a direct mode of differentiation using complex 
variables, as well as a previous implementation using handcoded direct differentiation for the mesh terms. Timings 
for the major components have been shown for several large-scale configurations, with dramatic savings obtained 
through the use of the adjoint formulation for the mesh linearizations. Demonstration optimizations yield improve- 
ments in vehicle performance for the discrete models used. 
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